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The wavelet transform, a family of orthonormal bases, is introduced as a technique for performing 
multiresolution analysis in statistical mechanics. The wavelet transform is a hierarchical technique 
designed to separate data sets into sets representing local averages and local differences. Although 
one-to-one transformations of data sets are possible, the advantage of the wavelet transform is as an 
approximation scheme for the efficient calculation of thermodynamic and ensemble properties. Even 
under the most drastic of approximations, the resulting errors in the values obtained for average 
absolute magnetization, free energy, and heat capacity are on the order of 10%, with a corresponding 
computational efficiency gain of two orders of magnitude for a system such as a 4 x 4 Ising lattice. In 
addition, the errors in the results tend toward zero in the neighborhood of fixed points, as determined 
by renormalization group theory. 

I. INTRODUCTION 

Spin models are popular tools for theoretical calculations and for numerical simulations, as their universality classes 
allow a huge range of different systems — as varied as binary metal alloys, surface adsorption, and neural networks — to 
be modeled simultaneously. For example, even the "trivial" one-dimensional Ising model can be used to model the 
helix-coil transition in biopolymers; the deep connection between magnetic models and polymer chains allow us to 
predict scaling behavior and other properties across an even wider range of materials.^ Lattice models are still widely 
used in modeling the thermodynamics of complex systems, because their regular structure simplifies the type and 
nature of interactions among components of the system. Moreover, the difficulty in obtaining analytical solutions 
of lattice systems, and the relative ease of computational simulations thereof, make them ideal test cases for new 
simulation algorithms. 

Although simulations of lattice models are relatively straightforward to implement, they share the same drawbacks 
as ofF-lattice models. The chief drawback is that as the number of particles grows large, the time required to sample 
the system accurately increases rapidly. A popular approach for addressing this problem is to coarse-grain the system: 
that is, we "rescale" the problem by increasing the basic size of a simulation element. For example, we might coarse- 
grain an atomic representation of a polymer chain into a "united atom" model, where a chain molecule is treated as 
if it consisted only of the backbone. More creative approaches redefine the problem to be addressed: for example, 
Mattice and coworkers have produced a method which maps a polymer chain atomistically onto a high-coordination 
lattice; this lattice is then used as the basis for a Monte Carlo simulation; the resulting configuration is then used to 
map back to continuous space to provide "evolution in time."^"^ 

This paper illustrates the use of the wavelet transform as a mathematical basis for performing thermodynamic 
computations of lattice models. The wavelet transform is an important tool in multiresolution analysis, which analyzes 
a system simultaneously at several length or frequency scales selected to reflect the actual physical processes underlying 
the observed behavior as closely as possible. The wavelet transform possesses a number of convenient properties, 
including orthogonality, compactness, and reconstruction; we will make these concepts more precise in Section II B 
below. The orthogonal nature of most wavelet constructions makes them a logical choice for use in ah initio density- 
matrix quantum chemistry computations, in which the selection of an accurate basis set is crucial to the convergence 
and efficiency of the calculations. Wavelet decompositions have been applied principally in electrical engineering, 
particularly in the field of signal processing. In this context, white noise and Markov processes have been studied using 
multiscale methods.^' To date, however, wavelet analysis does not seem to have been extensively applied to models 
in statistical mechanics. Huang uses wavelet analysis to observe the statistical distribution of multiplicity fiuctuations 
in a lattice gas,^^ while Gamero et al. employ wavelets to introduce their notion of multiresolution entropy, although 
their primary goal is dynamic signal analysis rather than statistical mechanics simulations,^^ while O'CarroU attempts 
to establish a theoretical foundation connecting wavelets to the block renormalization group. A more in-depth 
review of the connection between wavelets and renormalization theory is provided in a recent monograph by Battle. 
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II. WAVELET TRANSFORM FUNDAMENTALS 



A. The conceptual picture 

The wavelet transform is a hierarchical method for decomposing a data set into averages and differences. Like the 

Fourier transform, it can be used to provide a decomposition in both real space and reciprocal space (k-spacc), or time 

space and frequency space. Unlike the Fourier transform, however, it is capable of providing simultaneously localized 

transformations in both real and reciprocal space. A function localized in position space, such as a finite impulse 

function, cannot be represented by a few terms of its Fourier series: many terms are required before good convergence 

is achieved. By contrast, in wavelet space, this same function can be almost completely described by just a handful 

of wavelet coefficients. Although the first wavelet was discovered almost a century ago by Haar,^^ they have become 

an important computational technique only in the last decade, following the work of Mallat,^''' Daubechies,^^ and 
others. 20-22 

The wavelet transform, like any other transform, takes a mathematical object and transforms it into another: we 
can represent its action by writing 

w = W [«] ; (1) 

the specific form of W depends both on the type of wavelet we have selected, and the object u which we wish to 
transform. All versions of the wavelet transform W, however, are derived from the same source: a set of coefficients 
which define the transform. If u is a discrete data set, such as a signal sampled at regular intervals, then W is usually 
represented as a matrix; while if u represents a continuous data set, such as the same signal measured at all times, 
then yV acts as an integral operator. While the matrix form of W is often called a "filter bank" and the integral form 
a "wavelet transform," we will not distinguish between them in what follows, as the theory developed here for discrete 
lattices and filter banks should carry over to continuous systems and wavelet transforms essentially unchanged. 

Similar to the Fourier transform, the wavelet transform decomposes the object x into two separate components, as 
two different functions, a scaling function (j> and a wavelet function ip^ both operate on x. However, the two functions 
separate its components not into cosines and sines, but into averages and differences, with a "wavelength" equal to 
the "window" over which the scaling and wavelet functions are nonzero. In another important distinction, the wavelet 
transform is recursive^ so that it can be applied in succession to any set of averages which is produced using that 
wavelet transform, to produce another level of averages and another level of details. 

We shall now make the above concepts more mathematically precise. Let us define u to be a discrete set of samples 
u — (m (1) , M (2) , . . . , u (n)). Then applying the scaling and wavelet functions ^ and to u create a set of averages 
s {i) and a set of differences 5 (i): 

r-1 

s{i)=^cPik)u{i + k), (2) 

fc=0 

r-1 

6{i)=Y,i'ik)u{i + k), (3) 

fe=0 

where r is a finite integer which defines the length scale, often referred to as the "size of the support," over which 
(p and ijj are nonzero. The index i runs from 1 to n; generally the data set is padded with zeros to ensure that all 
sums in (2) and (3) are well-defined, although periodicity is sometimes used instead. ^'^ The coefficients (/>(fc) and 
ip (fc) in (2) and (3) are related, ^^'^^ and are central in controlling the features of the wavelet transform. Note the 
wavelet transform is inherently redundant: for every sample u (i) in the original set u, we now have two values, a 
local average s(i) and a local difference d{i). Since the new data are simply linear combinations of the original 
values, it is superfluous to retain both sets; at the same time, it is obvious that we cannot simply discard one set 
of data and recover all the original information using only the other data set. Instead, we choose to keep only 
the odd-numbered s(i)'s and (5(i)'s, eliminating the even-numbered samples; this process is called downs amplingP'^ 
Downsampling removes half of the s (i) and half of the S (i), regardless of the length r of the wavelet. Now we are 
left with n data: s (1) , s (3) , . . . , s (n — 1) and (5 (1) , (5 (3) , . . . , i5 (n — 1). These n data points can be stored as the 
level-one wavelet transform u oi u, by assigning s (1) , s (3) , . . . , s (n - 1) to {((1) (1) , {((1) (2) , . . . , (n/2), and the 
corresponding S (i)'s as u^-^'^ {n/2 -|- 1) , . . . , ii'^^^ (n). [The superscript (1) denotes that the wavelet transform has been 
applied once to this data set.] We can either stop at this level of description, or continue by further decomposing the 
averages: then the new object u^^^ to be transformed is u^^^ (1) = u (1) , . . . ,u^^^ (n/^) = {((n/2), and so on. Note 
that although u^^^ contains the averages s(i)'s and the differences ^(i)'s obtained in the previous step, successive 
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transforms only apply to averages obtained in the previous step. This process can be repeated until we have reduced 
our set of averages to a single point; no further averaging is possible. We assume henceforth that the data set u^''^ 
has been sufficiently downsampled to retain only the minimum data set required. 

B. Properties and examples of wavelet and scaling functions 

Until this point, we have not introduced any specific wavelet or scaling functions. Before we do so, we note that the 
choice of a wavelet transform to apply to a given system usually hinges on the desired properties which one wishes 
to include in the transformed data. Three principal properties are almost universally required for filter banks and 
wavelet families:^^'^^ 

1. Perfect reconstruction: No data is distorted by performing analysis followed by synthesis, so that the only 
permissible change is a delay in recovery of the original sample. 

2. Orthogonality: Wavelets computed at different length scales or at different spatial locations are mutually or- 
thogonal; thus fluctuations in the system are localized at the scales where they are most relevant. 

3. Compact support: Properly designed wavelets are identically zero except for a finite interval, which means that 
exact results can be obtained using only a finite number of terms. 

Other properties, such as orthonormality, symmetry in the functional form of the wavelet or a certain number of 
vanishing moments, can be taken into account when constructing the wavelet transform.^'* 

The two most commonly encountered selections arc the Haar and Daubechies wavelets, named after their respective 
discoverers. The Haar pair is the oldest and simplest set of wavelets:^^ the coefficients of the scaling function are 
(j) = {(j){0) (1)) = ^ (1, 1), while the coefficients of the wavelet function are ip = {ijj (0) , (1)) = (— 1, !)• No 
other wavelet can be described with two points, and therefore no other wavelet has a support as compact as the Haar 
wavelet. The scaling function (p simply averages the values stored at neighboring points, while tp finds the difference 
between those values; the extra factor of \/2 is incorporated to ensure orthonormality between overlapping 4> (k) and 
Ip (fc). A simple example of the action of the Haar wavelet is shown in Fig. 1. 

The Daubechies wavelets are a family of orthonormal functions whose construction was explicitly designed to 
have orthogonality as well as vanishing higher-order moments. Daubechies was able to show that the Haar wavelet 
is in fact the "first" member of the Daubechies family; that is, the Haar wavelet is the Daubechies wavelet with 
the shortest support. The second such member has four terms in its definition: the scaling function is defined by 
(j)= {(j){0) (1) ,(j){2),(j) (3)) = (1 + \/3, 3 + \/3, 3 - \/3, 1 - \/3) . The wavelet function reverses the order of the 
coefficients and inverts the sign of every other component, which allows the orthonormality properties to be satisfied: 
tjj = {ip{0) (1) ,ijj{2) ,tp (3)) = {—(p (3) ,(p{2) , —(p (1) , (p (0)). We can see that the Haar wavelet obeys the same 
pattern as the Daubechies wavelet: ipH = {—(p{l) ,(p{0)). This pattern can be extended, using different coefficients 
but the same general sign rules, for wavelets with 6, 8, 10, . . . coefficients. The resulting wavelet and scaling functions 
become increasingly smooth, and therefore are better suited for data sets in which there is only a gradual change in 
the data set with position — in thermodynamic systems, this would be more useful for, say, a spin- A'' Ising model than 
a spin- 5 Ising model (presuming that N ^ 

C. Matrix formulation of the wavelet transform 

For discrete systems, a conceptually simple method of implementing the wavelet transform is to set up the transform 
as a matrix equation. The input u is converted into a column vector u, so that the coefficients s (i) and 6 (?) are 
obtained via the dot product of u with vectors h (i) and 1 {i), where the vectors are padded so that the first nonzero 
element is located at position i: 

h(i) = (0,... ,0,<^(0), </)(!),... ,0(r-l),O,... ,0), 
l(i) = (0,... ,0,V(0),^(1),... ,V(r--l),0,... ,0). 

The vectors s and S, which contain the wavelet-transformed coefficients of the decomposition, can be obtained by 
forming the matrices H and L and right-multiplying by the vector u: 

s = Hu and 5 = Lu, 
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where the rows of the matrices H and L are the vectors (h(l),... ,h(n)) and (1(1),... ,l(n)), respectively. 
As mentioned in the previous section, we do not need to keep aU of the s(i)'s and i5(i)'s in order to ob- 
tain perfect reconstruction of our signal; thus, we can obtain all the necessary coefficients in a single ma- 
trix multiplication by combining the relevant rows of H and L into a single matrix W^^-*, whose rows are 
(h (1) , h (3) , . . . , h (n — 1) , 1 (1) , 1 (3) , . . . , 1 (n — 1)). Thus, the wavelet transformation (1) can be written as 



s(l) 




h(l) - 


s(n- 1) 




h(n- 1) 


S{1) 




1(1) 


(5(n- 1) 







nil) 
u{2) 



u (n) 



(4) 



We will denote the product on the left-hand side of (4) as u. 

As stated above, the wavelet process can be applied recursively: the set of averages (s (1) , s (3) , . . . , s (n — 1)) can be 
treated as a new data sample u^^^ = {u^^^ (1) , • • • , "^^^ ('^/2)), and operated on by an ^ x y reduction of W'^^, which 
we denote W^^^ to produce a new set of n/A averages (s*^^^ (1) , s^^^ (3) , . . . , s^^^ {n/2 — 1)) and corresponding new 
set of n/4 differences {5^^^ (1) , 5^"^^ (3) , . . . , J^^) (n/2 - 1)) . To reconstruct the original data set, we combine these n/2 
values along with the n/2 differences (^^^^ (1) , ^^^^ (3) , . . . , 5^^^ [n — 1)) obtained from applying W^^^. This process 
can be repeated as many times as desired, dividing the m non-downsamplcd averages s^*^) into m/2 averages s^'^"'"-'^) 
and m/2 differences S^'^^^h However, since at each iteration the matrix W^''^ only operates on selected elements of 
the vector u^*'^ = W^*^~^^u('^~^^, computations for multiple levels can be performed at the same time. Thus, if we 
wish to apply the wavelet transform K times, we can write this operation as an extended matrix product:^^ 



K 



uW=Wu=[]qWu, 
where the Q^*^' are a family of matrices of the form 



w^'^) 

I 



(5) 



(6) 



In (6), Q*^*^^ is always a,n N x N matrix, while the matrix W^*) has size {N/2 



k-l\ 



{N/2 



To recover the original data sample u following a wavelet transform W, we can simply multiply u by the inverse of the 
wavelet matrix W. The matrix W is unitary: that is, its inverse is equal to its transpose W-^. Consequently, 

if W is known, all that is necessary to reverse the transformation is to left-multiply u'^'^' by the transpose W'^. 
Moreover, the computation (4) of the wavelet transform can usually be carried out "in place" by manipulating local 
coordinates; in this manner, the computation is carried out even more rapidly than a standard multiplication, and 
without the increased storage costs associated with matrix multiplications.^^' 



D. Multidimensional wavelet transforms 



Since virtually all problems in lattice thermodynamics are in multiple dimensions, it is necessary to take the 
wavelet transform of a multivariate function or data set. Several methods have been developed to carry out such 
transformations; among them are Cohen and Daubechies's separable wavelets, which form the multidimensional scaling 
and wavelet functions (j){x,y), tj)xx{x,y), ipxy{x,y), and ^pyy{x,y) from products of the one-dimensional scaling 
and wavelet fimctions (l){x) and xp(x).''^'^ A more general algorithm, the lifting algorithm,, has been developed by 
Sweldens.^^'^° It divides the wavelet transform into two steps: the first computes the wavelet coefficients 6 {i) ; the 
second step uses the wavelet coefficients to speed up the calculation of the scaling coefficients. "Initialization" of the 
lifting algorithm requires the use of an appropriately selected basis function. 

A particularly convenient basis function for the multidimensional lifting transform is the generalized orthogonal 
Haar wavelets outlined by Sweldens.^® An extension of the one-dimensional wavelet transform, they can be cre- 
ated in any number of dimensions, and have the same basic orthonormality properties as the one-dimensional Haar 
functions, although the orthonormality constant becomes 2""^/^, where d is the dimensionality of the system. [The 
two-dimensional version is shown in Figure 2.] Moreover, the use of the Haar wavelets as a starting point for further 
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iterations of the lifting algorithm allow the development of additional, "better" lifted wavelets with more desirable 
properties, such as smoothness. 

It can further be seen that the "oversampling" problem which exists in one dimension will be magnified in multiple 
dimensions: since the number of wavelet functions which are produced from each data point increases by a factor of 
two with each additional dimension added, we must reduce the number of points maintained for each wavelet function 
by that same factor. Thus, in two dimensions, we keep only every fourth point; in three dimensions, every eighth 
point, and so on. The wavelet transform for multidimensional systems can thus still be written in the form of (5) and 
(6), after we have written the multidimensional data set in terms of a column vector. This can be accomplished by 
wrapping around the edges of the matrix in creating u: for example, after inserting element (1, N) of a two-dimensional 
data set into u, we next store element (2, 1), and so forth. The other significant difference in the structure of these 
equations is that the size of the submatrix W^'') in (6) is now N /2^^-'^'^'^>(. instead of x 

III. WAVELET ANALYSIS OF LATTICE THERMODYNAMICS 

A. Applying the wavelet transform to Hamiltonians 

The standard model for studying the thermodynamic behavior of lattice systems is the spin-^ Ising model, which 
contains both nearest-neighbor pairwise interactions as well as interactions between lattice sites and an external field. 
The Hamiltonian for this system is normally written in the form 

-l3n = h'Y^ai + J^a.Gj, (7) 

i {ij) 

where h is the strength of the external field in the direction of the spins crj, and J is the strength of the interaction 
between nearest-neighbor pairs of spins on the lattice; these pairs are indicated by the subscript (ij) in the second 
summation in (7). The inverse temperature /3 = (/cbT)"^; for convenience we let ks = 1, so that temperature, 
external field, and nearest-neighbor interactions are all dimensionless quantities. The model can be further extended 
by the inclusion of a position-dependent external field, or by the inclusion of pairwise interactions beyond nearest 
neighbors; in this case, (7) can be written in the more general form 

-/3n = ^ /ijO-j + X! X] Jij^i^^J' (8) 

i i j 

where in (8) hi is the strength of the external field at lattice site i, and is the strength of the interaction between 
sites i and j. Analytical solutions of (7) are well-known for one-dimensional systems and two-dimensional systems 
when h — 0;'^^'"'^ it has recently been suggested that analytic solutions do not exist for more complicated systems.''"^ 
While (7) and (8) are compact representations of the Hamiltonian of the system, the expansion of the lattice 
variables Si and Sj as a sum of wavelet coefficients makes these equations impractical for applying the wavelet 
transformation. Since the system is described discretely, we want to use discrete wavelets, and therefore a matrix 
formulation of the Hamiltonian would be convenient. Using graph theory,'^'* this is readily accomplished: let the 
vectors u = (cti, (T2, - ■ ■ , <7n) and h = {hi, . . . , /iat) denote the values of each of the N lattice variables in the system 
and the set of external-field interaction strengths, respectively, constructed in the row-wise manner described in the 
previous section. Furthermore, define the matrix J such that element Jij is the strength of the interaction between 
site i and site j. If these sites do not interact, then Jij = 0. Then, the Hamiltonian (8) can be written in the form of 
a matrix equation: 

-/m = h^u + u^Ju, (9) 

where the superscript T denotes the transpose of the vector (or matrix) which precedes it. 

The matrix W which defines the wavelet transform satisfies by construction W-^W = I, where I is the identity 
matrix. Therefore, to apply the wavelet transform, we simply insert W-^W between each pair of terms in (9), thereby 
obtaining 

-PH = (h^W^) (Wu) + (u^'W^') (WJW^) (Wu) . (10) 



5 



Using the general matrix property that B-^A-^ = (AB) , (10) can be rewritten in terms of the wavelet-transformed 
vectors h(^)= Wh, u(-^)= Wu, and the wavelet-transformed matrix 5^-^)= WJW"^: 

-m = (h(^))^uW + (uW)^ jWu(^). (11) 

It is important to note that in writing equation (11), we have not made any expUcit assumptions about the form of the 
matrix W, other than to require it to be a matrix describing a wavelet transform. As a result, wc can perform several 
levels of multiresolution simultaneously through appropriate preparation of the matrix W in the manner outlined 
above. Using the inverse wavelet transform = to recover an original configuration u is a one-to-one 

mapping only if wc arc provided with n distinct data points, as contained in u^^^. 

Examining (11), we see that the result of applying the wavelet transform to a set of spins u is to create a new 
representation u^^) , which contains the same information about the spins as does the original state vector u. However, 
the vector u^^' contains averages, where d is the dimensionality of the lattice; these averages can be viewed as 

"block spins" in a sense similar to that of Kadanoff."^^ The remaining elements of u*^''^^ contain the local differences in 
the spins; that is, they can be used to describe the specific set of spins which give rise to a particular block spin s-^\ 



B. Computing thermodynamic functions 

The canonical partition function 

Z = ^exp(-/3W(u)), (12) 

where S is the configuration space of the system, can be used to derive all the thermodynamic properties of a lattice 
system. Applying the wavelet transform W to the state vector u results in a new state vector u^^^ belonging to the 
configuration space S^^-*. Provided that W satisfies the perfect reconstruction property, if the summation over u G S 
in (12) is replaced with the summation over u^-'^) g S*^''^^ the results will be identical. This result follows naturally, 
since perfect reconstruction necessarily implies that there is a unique state vector u^^^ e §^^^ for each state vector 
u G §, and by construction of the wavelet-transformed Hamiltonian (11), = — /JTi^^'. 

Ensemble averages involving wavelet-transformed variables are in general no more complicated than computations 
involving the original variables. In general, the transformation of a function (or functional) of the original variables 
/ (u) or / [u] will generally have the same characteristics after applying the wavelet transform to obtain f^^^ (u*^-^)) 

or /(^) [w^^)] . Moreover, the standard properties of ensemble averaging, such as linearity, also apply, which makes 
calculations of moments of the distribution particularly simple. 

The above formulation can be applied to any system whose Hamiltonian can be written in the form of equation 
(8), or as the sum of contributions, each of which is of that form. While the wavelet transform can be applied to 
any lattice system whose Hamiltonian is a function of the components of the state vector u — or, indeed, to any 
Hamiltonian which is a functional of u (r) for continuous systems — in most of these cases, it is necessary to rely on 
the more cumbersome series expansions, or integral transforms in the case of continuous systems. 

In addition, for Ising spin variables, the use of the wavelet transform poses an additional challenge. While it is 
entirely straightforward to describe the possible values an individiial lattice site can take — for a spin-g Ising model, the 
allowed values of an individual spin a are —q,~-q+l,...,q — l,q —the rules which determine whether an arbitrarily 
selected "transformed" state vector u*^-^) represents a real state vector u are cumbersome to manipulate and to 
implement, and have proven a formidable challenge in prior research as well.^^ 



IV. ANALYSIS 



In this section, wc use two variants of the two-dimensional Ising model as a basis for our calculations: wc look at 
4x4 and 32 x 32 Ising lattices; the former is used for calculations when it is desired to perform calculations over 
all states explicitly, while the latter is illustrative of larger systems, for which exact calculations are intractable. The 
emphasis in this paper will be on the use of wavelets to yield approximate answers in significantly faster time than 
is possible with a Monte Carlo simulation incorporating all degrees of freedom explicitly. The methodology by which 
the wavelet transform can be extended to Monte Carlo simulations of lattice systems is the focus of the companion 
paper. 
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Before proceeding to the results of the calculations, we make note of the time required to execute the simulations. 
Each of Figures 3 through 6 plots the observed variables in the temperature range T = 0.50 to T = 5.00, with a step 
size of AT = 0.01. The computations for the original problem, with 2^® = 65536 states to consider explicitly, required 
more than 6.8 seconds per point to execute the required calculations on a 733 MHz Pentium III; the same calculations 
using one and two applications of the wavelet transform required just 0.061 and 0.026 seconds per point to consider 
5* = 625 and 17 states, respectively. 

A. Weighting functions for wavelet-transformed statistics 

Until now, we have worked with wavelet transforms which preserve the number of degrees of freedom between 
the original and transformed problems. This approach yields results for the wavelet-transformed system in exact 
agreement with those for the original system. However, as mentioned above, in such cases the transformed equations 
are usually harder to model than the original ones. Hence it is desirable to use the wavelet transform not only as a 
means of describing a lattice system, but also to derive an approximation scheme whereby estimates of thermodynamic 
properties can be made efficiently, while still offering error estimates that bound the true results. 

The approach we adopt in this paper is to ignore all local differences: that is, we assume that S^''^ = for all values 
of i and k. As a secondary assumption, we assume that correlation functions which include wavelet coefficients are 

also equal to zero: that is, we assume (^S^'^'' A {■)^ = for any choice of property A{-). These extreme assumptions 

represent a "worst-case scenario" for the use of the wavelet transform method; any more accurate representation of 
the behavior of the wavelet coefficients will lead to similarly more accurate results in the calculations we present 
below. Note that under certain circumstances, these approximations are accurate: for Hamiltonians of the form (8) 
which do not exhibit interactions with an external field, there exist configurations with equal energies but opposite 
signs for 5- ; consequently, when we take the average over all configurations, the ensemble average of will vanish. 

Let us assume that we have applied the wavelet transform to our original Hamiltonian, Ti. (s) , and have obtained a 
new function (flW). If we then make our approximation that all the wavelet coefficients are negligible, then we 

can reduce H^^^ (u^-'^^) to a new function H'^^^ (s'-''^^-'), where s^-*^^ represents the set of all averages (^s'{^\ . . . , s^^^ 

that were preserved by the wavelet transform. We know that the partition function of the system before we performed 
the wavelet transform is given by (12); after using the wavelet transform, we hypothesize that the partition function 

of the new system is given by 

Z= Yl (s^""^) exp (s(^))) , (13) 

s(K)gs(/f) 

where the w (s^^^) is the weight of configuration s*^^) in the configuration space S^^'; since multiple configurations 

u G S can correspond to the same s*^^\ we cannot set w (s^^^) = 1 as we did in the untransformed Ising model. 
Clearly, the new partition function (13) will be identical to the original partition function (12) if we define 

w{s)= Yl exp(-/3(H(u)-HW(sW))), (14) 

u:u— >s(^) 

where the notation u : u ^ s^^^ denotes that the sum is to be performed over all configurations u which have the 
same set of wavelet-transformed averages s^''^-'. However, evaluating (14) to obtain the weighting functions is no more 
tractable than the computation of the original partition function. Thus, any computational efficiency to be gained is 
by finding an economical approximation for (14). One such approach is to take w (s^^^) to be equal to the number 

of states s whose wavelet transform yields s^^^, so that (13) is the standard form of the canonical partition function 
for a system with degenerate energy levels. Let us call the number of such states s with equivalent averages s^^^ the 
degeneracy of state s^^\ and denote this degeneracy as g (s'^') . The restriction of u^^' to the averages s^^^ prevents 
a unique reconstruction of u, unless g (s^^^) = 1. 
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B. Order parameter 



A natural variable to compute is the order parameter, generically denoted r], which for lattice spin systems is the 
average magnetization, 




for other members of the Ising universality class, the order parameter can represent the overall density p or the 
difference between the densities of two phases.'^® While the computation of the order parameter is straightforward 
in simulations, its calculation can be made difficult in the case of zero external field because of symmetries in the 
configuration space: for every configuration with magnetization m^, there exists another configuration with the same 
energy and magnetization —rrii. As a result, when all the states are combined using (15), we do not observe spontaneous 
magnetization, but instead find m = at all temperatures. Thus, we consider the absolute value of the magnetization, 
|m| in place of the magnetization m. The results of this calculation for the 4x4 Ising model with the wavelet transform, 
setting all S (i)'s to zero, and without the wavelet transform are shown in Figure 3. We note that the error is essentially 
negligible for temperatures below T = 1, and decreases again for large values of T, where differences in energy levels 
become negligible and the average magnetization of a state is the only contributing factor to |m|. 



C. Free-energy considerations 

The Hclmholtz free energy is just the logarithm of the partition function: A = —ksThiZ] if we estimate the 
partition function using an expression such as (13), wo naturally expect the approximated value A to differ from 
its true value. When we examine the behavior of the 4x4 Ising model under the wavelet transform, we find as 
expected that the two free energy surfaces are similar, although they are clearly not identical. In particular, the exact 
numerical values obtained from the two equations are not the same. However, since the assignment ^ = is arbitrary 
in any system, we can choose to define A = as cither the maximum or the minimum free energy obtained in each 
system. Under these conditions, the energy scales for the exact calculation and the wavelet transform calculation 
are essentially identical, particularly at the so-called "fixed points" of the system — that is, for points which are not 
affected by rcnormalization transformations.^^ For the two-dimensional Ising model, these points are at T = 0, and 
at infinite (positive or negative) values of the external field interaction h. 

This agreement at the fixed points should not be surprising: the fixed points correspond to unique configurations of 
the system, such that for the configurations s*^^) which result from the wavelet transform of the fixed points u* , the 
degeneracy of the states is unity. Consequently, the approximation w (s*-^-*) = g (s'^-*) is correct for the dominant 
configuration in the system, and therefore the behavior of the approximate partition function Z' is almost identical 
to that of the true partition function Z in the vicinity of the fixed points. 

As can be seen in Figure 4, the free energy converges to the same values in the low-temperature limit, where only a 
few states which are essentially unaffected by the wavelet transform contribute to the partition function of the system. 
At high temperatures, the agreement is less exact, because of the approximations made in evaluating the Boltzmann 
factors of the block spin configurations. 



D. Entropy of a wavelet-transformed lattice system 

It is relatively straightforward to show that coarse-graining a system with a countable phase space — that is, a phase 
space whose states can be completely enumerated — requires a correction factor to ensure satisfactory agreement with 
results from the fine-scale system. As an example, consider the entropy of a wavelet-transformed spin-^ Ising model. 
The original lattice has 2^' total configurations, and therefore in the high-temperature limit, the entropy approaches 

S^^^ = kBNt\n2. (16) 
If we ignore the weighting factor of the wavelet-transformed system by setting w (s*^^^) = 1 for all states i*^^', the 
resulting system will have ( 1 + ni=i ) ' the maximum possible entropy for this system is 

5i^ = fcBM^'ln(^n'^'''+l)- (17) 
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Comparing (17) to the original limit of S'max = kBNtln2, we must have Smax > S'max, since the configuration space is 
smaller, which indicates that our coarse-graining has effectively reduced the available entropy of the system. We can 
conclude from this that unless the volume of phase space is conserved by the coarse-graining procedure, some entropy 
will be lost at high temperatures, regardless of the accuracy of the coarse-graining. Our choice of the degeneracy 
g (s^^^^) for the weighting function w (s*-^^) preserves the volume of phase space in the limit of T ^ oo for the systems 
under study here. 

Looking at the specific case of the 4x4 Ising model, we compute the entropy as a function of temperature for 
zero external field for the exact problem, and for the wavelet-transformed problem using one and two iterations of 
the two-dimensional Haar wavelet. The results are shown in Figure 5. As before, we find that errors vanish in the 
low-temperature limit, where only the lowest-energy states make a contribution to the partition function; since there 
are two such states, namely those with all spins up and all spins down, we find that the T — > limit of the entropy is 

5(T-0)«-^ilni=ln2 

i=l 

In the high- temperature limit, the entropies agree because, by construction, we have X)ses^(^) ~ 2^^, and therefore, 
since the Boltzmann factor goes to unity for all states as T ^ oo, Z ^ X^gei^ (^)' entropy tends toward 

the value S'max given by (16). Disagreement in the intermediate temperature regime is largely the result of grouping 
together states with different energies into a single transformed state with a single energy. In particular, the grouping 
of higher-energy states together at a lower energy increases the probability of being in those states, and therefore 
increases the total entropy of the system. This explains the seemingly anomalous result of coarse-graining increasing 
the entropy at intermediate temperatures observed in Figure 5. 

For more general systems, we can consider the formal definition of the entropy for a continuous distribution, 

S = -kB [ dup{u)lnp{u) , (18) 

where p (u) is the probability of observing the configuration u. For the original Ising model, each configuration u is 
unique, and therefore the probability p (u) is just the Boltzmann weight exp {—(3TI (u)) , and thus the entropy is 

given by 

S = kBlnZ+^ [ rfu/3W(u)exp(-/3W(u)). (19) 

After performing the wavelet transform, and discarding finer-scale details, we must account for the weighting factor 
w (s) in our expression for the entropy. Consequently, the entropy for such a system is given by 

= - /cb / dsw{s)p{s)\np{s) 

= kBlnZ+^[ ds«;(§)/3WW(sW)e-^«'"H^'"'), (20) 
Z Jses ^ ' 

where the coarse-grained partition function Z is given by (13). 



E. Fluctuation properties 

For a given thermodynamic system, the constant-volume heat capacity C is defined as the fluctuation in the internal 

energy of the system: 

It is well-known that in the vicinity of a continuous phase transition, the heat capacity C diverges. At the same time^ 
it is also known that no system of finite size can have C ^ oo, since the energy, and hence the variance (-E^) — {E) 
in the energy of the system will also remain finite. However, we can still observe evidence of power-law divergence 
near the critical temperature.^" For small finite systems, we do not observe evidence of a divergence in the heat 
capacity; instead, the heat capacity is a smooth function of temperature T. The wavelet transform largely preserves 
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the behavior of the original model; we observe the same general functional form in both the original model as well as 
the transformed systems, as seen in Figure 6. However, for large finite systems, we can still detect the characteristic 
power-law divergence in the vicinity of the critical point; such a divergence is shown for a 32 x 32 Ising model as 
Figure 7. 

Applying the wavelet transform to a given thermodynamic system, we expect that we will still find evidence of a 
critical point; however, because of the mean-field like behavior of the wavelet transform method, the critical point 
should be found at a higher temperature in the wavelet-transformed system than in the original system. Computing 
the heat capacity for large lattices cannot be done by exhaustive enumeration; therefore, we save discussion of our 
numerical results for a subsequent paper,'''' and present a brief heuristic argument. As the size of a block spin 
increases, corresponding to a greater number of iterations of the wavelet transform, we expect that the suppression 
of fluctuations in the system will lead to an increase in the temperature of the phase transition, while the actual 
numerical value of the heat capacity itself decreases. The decrease in the number of degrees of freedom in the system 
is accompanied by an increasing trend towards homogeneity of the system — the remaining configurations begin to look 
more and more similar to one another. As a result, the variance measured by (i?^) — {E) decreases with increasing 
amounts of coarse-graining. This trend is much less pronounced for small systems, as shown in Figure 6. When 
the number of degrees of freedom is small, the ability to perform an exact enumeration ensures that the only errors 
observed are those introduced by the wavelet approximation itself. Thus, since the free energy and entropy of small 
systems are accurately reproduced by the transform, only minor deviations in the phase transition temperature for 
small systems are expected. 

V. SIMILARITIES TO RENORMALIZATION GROUP THEORIES 

It should already be apparent that close aSinities exist between the method described here and position-space 
renormalization methods. The construction of wavelet-based averages naturally corresponds to Kadanoff's concept 
of a "block spin" transformation:''^ both the present method and Kadanoff's approach rely on combining a region 
of contiguous spins into a single new "block spin." However, with the approach outlined above, we do not seek to 
impose the requirement that the block spins must be restricted to have the same spin values as the original spins. 
For example, in a wavelet-transformed spin-i Ising model, the block spins can take on values othca- than -1-1 and —1. 
In addition, we do not impose "majority rules" or other "tiebreakers" in the case where there arc equal numbers of 
"up" and "down" spins contained in a single block. 

A less apparent connection can also be established with Migdal's "bond-moving" approximation,''^' or similarly 
with Wilson's recursion method.*^ In particular, Migdal's method can be compared to applying a separable wavelet 
transformation, in that the recursion is performed in only one spatial direction at a time. However, using the methods 
of either Migdal or Wilson, after performing the bond-moving or decimation transformation, the new interaction 
strengths J-j are determined by manipulating the resulting Hamiltonian and recasting it in the form of the original, 
leading at times to very complicated, nonlinear formulas for the J^j as a function of Jij. By contrast, using wavelets, 
the values for the transformed coupling coefficients are obtained directly by transforming the interaction matrix J, 
since J = WJW^. The trade-off that must be made for the algorithmic transparency of obtaining the new coupling 
constants via a matrix multiplication is that we cannot solve for the fixed points of the recursion relation, which 
makes the production of a renormalization "flow diagram" using the wavelet method a difficult problem. 

The behavior of the wavelet transform and of the renormalization group differs in another, important manner. 
By the nature of its construction, the group of renormalization transformations is only a semi-group, inasmuch as 
reversing the mapping to move from a coarsened description to a more detailed description is impossible. Using 
the wavelet transform, a reverse mapping is theoretically possible: reversibility fails because of the approximations 
invoked to reduce the number of wavelet coefficients which are kept, rather than as an inherent limitation of the 
wavelet transform itself. 



VI. CONCLUSIONS 

Using the wavelet transform as a basis for calculations of lattice systems yields an impressive reduction in the 
computational time required to obtain estimates of thermodynamic properties, at the price of modest errors in 
accuracy, given the relatively small size of the systems considered here, and the approximations introduced to simplify 
our calculations. The wavelet transform provides a systematic approach to coarse-grain systems to any level of 
resolution. The method produces correct results in the vicinity of fixed attractors, with decreasing accuracy as one 
approaches the critical point in parameter space. 
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In addition, the use of the wavelet transform on such systems is simple to implement: after selecting a particular set of 
wavelets, all of the operations can be reduced to matrix multiplications, as shown in Section III A. The computational 
requirements needed to implement the transformation are small: if the transformation cannot be accomplished in the 
explicit form of a matrix multiplication, it is possible to organize the calculation to be performed in-place.^^ Moreover, 
the amount of coarse-graining that can be achieved using the wavelet transform can be adjusted dynamically. These 
properties of the wavelet transform lead naturally to an extension of the method to systems for which simulations are 
required; the resulting algorithm is the focus of the companion paper. 



VII. ACKNOWLEDGMENTS 



Funding for this research was provided in part by a Computational Sciences Graduate Fellowship (AEI) sponsored 
by the Krell Institute and the Department of Energy. 



^ K. F. Freed, Renormahzation Group Theory of Macrornolecules, Wiley, New York, 1987. 

^ R. F. Rapold and W. L. Mattice, Macrornolecules 29, 2457 (1996). 

^ J. Cho and W. L. Mattice, Macrornolecules 30, 637 (1997). 

" P. DORUKER and W. L. Mattice, Macrornolecules 30, 5520 (1997). 

^ P. DORUKER and W. L. Mattice, Macrornol. Theory Simul. 8, 463 (1999). 

'^T. A. Arias, Rev. Mod. Phys. 71, 267 (1999). 

G. Beylkin, N. Coult, and M. J. Mohlenkamp, J. Comput. Phys. 152, 32 (1999). 
* B. R. Johnson, J. P. Modisette, P.J. Nordlander, and J. L. Kinsey, J. Chem. Phys. 110, 8309 (1999). 
^M. R. Luettgen, W. C. Karl, A. S. Willsky, and R. R. Tenney, IEEE Trans. Sign. Process. 41, 3377 (1993). 

M. R. Luettgen, Image Processing with Multiscale Stochastic Models, PhD thesis, Massachusetts Institute of Technology, 

1993. 

" D.-W. Huang, Phys. Rev. D 56, 3961 (1997). 

L. G. Gamero, a. Plastino, and M. E. Torres, Physica A 246, 487 (1997). 
" M. O'Carroll, J. Stat. Phys. 73, 945 (1993). 
" M. O'Carroll, J. Stat. Phys. 71, 415 (1993). 

G. Battle, Wavelets and Renormalization, World Scientific, Singapore, 1999. 

A. Haar, Math. Ann. 69, 331 (1910). 
" S. G. Mallat, IEEE Trans. Pattern Analysis and Machine Intelligence 11, 674 (1989). 
1* S. G. Mallat, Trans. Amer. Math. Soc. 315, 69 (1989). 

I. Daubechies, Ten Lectures on Wavelets, volume 61 of CBMS-NSF Regional Conference Series in Applied Mathematics, 

SIAM, Philadelphia, 1992. 
^° A. Cohen and I. Daubechies, Rev. Mat. Ibero-amer. 9, 51 (1993). 

W. SWELDENS and R. Piessens, SIAM J. Numer. Anal. 31, 1240 (1994). 

W. SwELDENS and R. Piessens, Numer. Math. 68, 377 (1994). 

G. Strang and T. Nguyen, Wavelets and Filter Banks, Wellesley-Cambridge, Cambridge, MA, 1996. 
W. Sweldens and P. Schroder, Building your own wavelets at home, in Wavelets in Computer Graphics, pp. 15-87, ACM 
SIGGRAPH Course notes. New Orleans, 1996. 
"^^ I. Daubechies, Commun. Pure Appl. Math. 41, 909 (1988). 

D. GiNES, G. Beylkin, and J. Dunn, Appl. Comput. Harmon. A. 5, 156 (1998). 

G. H. GOLUB and C. F. Van Loan, Matrix Computations, Johns Hopkins University Press, Baltimore, 1996. 

L. N. Trefethen and D. Bau III, Numerical Linear Algebra, SIAM, Philadelphia, 1997. 
2^ W. Sweldens, SIAM J. Math. Anal. 29, 511 (1997). 

I. Daubechies and W. Sweldens, J. Fourier Anal. Appl. 4, 247 (1998). 
'^^ R. K. Pathria, Statistical Mechanics, Butterworth-Heinemann, Woburn, MA, 1996. 
^2 L. Onsager, Phys. Rev. 65, 117 (1944). 

S. ISTRAIL, Statistical Mechanics, Three-Dimensionality, and NP-Completeness: I. Universality of Intractability for the 

Partition Functionof the Ising Model Across Non-Planar Surfaces, in Proceedings of the 32nd Annual ACM Symposium on 

Theory of Computing, pp. 87-96, Portland, OR, 2000, ACM Press. 

T. H. Gormen, C. E. Leiserson, and R. L. Rivest, Introduction to Algorithms, McGraw Hill-MIT Press, Cambridge, 
MA, 1990. 



11 



N. GOLDENFELD, Lectures on Phase Transitions and the Renormalization Group., Addison- Wesley, Reading, MA, 1992. 

C. Best, A. Schafer, and W. Greiner, Nucl. Phys. B: Proc. Suppl. 34, 780 (1994). 

A. E. Ismail, G. Stepiianopoulos, and G. C. Rutledge, Submitted to J. Chem. Phys. (2002). 

H. E. Stanley, Introduction to Phase Transitions and Critical Phenomena, Clarendon Press-Oxford, Oxford, 1971. 

S.-K. Ma, Modem Theory of Critical Phenomena, Benjamin/Cummings, Reading, MA, 1976. 

J. J. Binney, N. J. DowRiCK, A. J. Fisher, and M. E. J. Newman, The Theory of Critical Phenomena: An Introduction 

to the Renormalization Group, Oxford University Press, Oxford, 1993. 

A. A. MiGDAL, Sou. Phys. -JETP 42, 413 (1976). 

A. A. MiGDAL, Sov. Phys.-JETP 42, 743 (1976). 

K. J. Wilson and J. Kogut, Phys. Rep. 12, 75 (1974). 

P.-G. DE Gennes, Scaling Concepts m Polymer Physics, Cornell University Press, New York, 1979. 
M. E. Fisher, Rev. Mod. Phys. 70, 653 (1998). 



12 



FIG. 1. (a) A sample signal u (n). (b, c) The one-dimensional Haar scaling function 4> (x) and wavelet (differencing) function 
i}{x). (d, e) The first-level scaling coefficients s{n) and 5{n) produced from the original signal u(n) using the Haar pair, 
following downsampling of the signal. 

FIG. 2. The two-dimensional orthogonal Haar wavelets. The coeSicient in all cases is 1/2, times the sign indicated in each 
quadrant. Note that the scaling function <f> is built from the one-dimensional Haar scaling function, while the wavelet functions 
are built from the one-dimensional wavelet function. 

FIG. 3. (a) Average absolute magnetization |m| = \ai\ of the 4x4 Ising model at zero external field as a function of 
temperature, as computed by an exact enumeration using no wavelet transform (solid line), and using one and two iterations of 

the two-dimensional Haar wavelet (dot-dashed and dotted lines, respectively) . (b) Error in the average absolute magnetization 
of the 4x4 Ising model for one and two iterations of the two-dimensional Haar wavelet versus an exact enumeration. 

FIG. 4. Free energy A of the 4x4 Ising model at zero external field as a function of temperature, as computed by an exact 
enumeration using no wavelet transform (solid line), and using one and two iterations of the two-dimensional Haar wavelet 
(dashed and dot-dashed lines, respectively). 

FIG. 5. Entropy of the 4x4 Ising model at zero external field as a function of temperature, as computed by an exact 
enumeration using no wavelet transform (solid line), and using one and two iterations of the two-dimensional Haar wavelet 
(dashed and dot-dashed lines, respectively). The bottom of the y-axis corresponds to the zero-temperature limit of S = In 2. 

FIG. 6. Heat capacity of the 4x4 Ising model at zero external field as a function of temperature, as computed by an exact 
enumeration using no wavelet transform (solid line), and using one and two iterations of the two-dimensional Haar wavelet 
(dashed and dot-dashed lines, respectively). 

FIG. 7. Heat capacity of a 32 x 32 Ising model, illustrating the power-law divergence in the vicinity of the critical point 
(identified by the arrow). 
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